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ABSTRACT 

We study the dynamics of the chiral phase transition expected during the expansion 
of the quark-gluon plasma produced in a high energy hadron or heavy ion collision 
in the 0(4) linear sigma model to leading order in a large N expansion for strong 
coupling constants. Starting from an approximate equilibrium configuration at an 
initial proper time T in the disordered phase, we study the transition to the ordered 
broken symmetry phase as the system expands and cools. We give results for the 
proper time evolution of the effective pion mass, the order parameter < (7 > as 
well as for the pion two point correlation function expressed in terms of a time 
dependent phase space number density and pair correlation density. We investigate 
the possibility of disoriented chiral condensate being formed during the expansion. 
In order to create large domains of disoriented chiral condensates low-momentum 
instabilities have to last for long enough periods of time. Our simulations show 
that instabilities that are formed during the initial stages of the expansion survive 
for proper times that are at most 3 fm/c. 



1. Introduction 

The possibility of producing large correlated regions of quark condensate < qtqj > 
pointing along the wrong direction in isospin space was proposed 1 to explain rare 
events where there is a deficit or excess of neutral pions observed in cosmic ray exper- 
iments 4 . The idea that such disoriented chiral condensates (DCC's) can be formed 
has been the source of several experimental proposals 2 ' 5 in high energy collisions and 
a number of theoretical papers 3 ' 6 ~ 12 . Theoretically, it has been recognized that the 
possibility of forming large regions of DCC relies on the existence of a substantially 
large regime in which the hot plasma formed after the collision evolves out of equilib- 
rium 3 . In fact, if thermal equilibrium is approximately preserved by the dynamics, 
the typical correlation length would be determined by the pion mass and therefore 
would be too small to matter. For this reason, a number of authors made several 
attempts to analyze the nonequilibrium aspects of the dynamics of the chiral phase 
transition. These approaches vary in form and content: some authors performed nu- 
merical simulations on classical models 3 > 9 > 10 5 others used phenomenological terms - 



inspired by the classical kinetic theory - to model the interaction between the con- 
densate and the quasiparticles 11 while some attempted to incorporate quantum and 
thermal fluctuations 6 ' 12 into the picture. 

There are clearly two important questions that the theoretical models should aim 
to answer. The first one is to determine whether during the evolution that follows 
the collision there are instabilities affecting the fluctuations, in which case there is a 
chance of the correlations growing. If this occurs, structure may tend to form through 
a process like spinodal decomposition. 

The second question is, assuming the instability exists, can the correlated domains 
grow large enough so that many pions can be emitted from each domain making the 
detection of DCC's possible. The time scale for the instabilities to exist is strongly 
influenced by the strength of the interactions. In the strong coupling regime insta- 
bilities exist only for a short period of proper time. This makes it difficult to obtain 
large sizes for the correlated domains. On the other hand in a weakly coupled system 
the domains can grow for long periods of time, and a different picture emerges. 

In our approach we do not put in instabilities by hand but look at typical fluc- 
tuations of initial conditions starting in a region of stability. Previous researchers 
have assumed the existence of an initial instability and focused their calculations on 
studying the growth of correlations. The typical and quite drastic way in which the 
instability had been previously introduced in the problem was by using the so called 
quench approximation. In a quench, the system (the hot plasma created after the 
collision) is subject to a sudden external action (the expansion into the vacuum) that 
has the following net effect: it does not change the state of the system but, instan- 
taneously modifies the effective potential turning it "upside down". This is clearly 
an idealization that may be entirely inappropriate. In this paper we will present an 
approach that does not require this approximation and that enables us to study (in 
a simplified model) the existence of instabilities in a rather straightforward way. 

Thus, our aim is to study the evolution of DCC's without imposing the quench 
approximation or any other ad hoc way of modeling the appearance of an instability. 
We start at an initial value of the proper time r in the symmetric phase, where the 
particle masses are positive, by choosing a thermal distribution of particles above the 
critical temperature T 2 = 12 f%/N 20 . The system then cools by expanding into the 
vacuum. The expansion can be described by the Bjorken boost invariant picture for 
cooling 13 that has been applied in various hydrodynamic models to hadronic as well 
as heavy ion collisions. 

The natural expansion and subsequent cooling of the plasma into the ordered 
vacuum is studied numerically by solving the update equations for the quantum modes 
as well as for the expectation values evolving in proper time. In this way, we are able 
to study the evolution of the plasma in a self consistent way, without imposing any 
instability by hand. We analyze various reasonable initial conditions on the fields, and 
determine whether they lead to instabilities. When the effective pion mass becomes 
negative instabilities ensue. We also determine the time evolving order parameter 
< a > and the spatial correlation function for the pion field. 

In the investigations done so far, simple phenomenological models have been used 



hoping that they describe the fundamental physics involved in the dynamics of the 
phase transition. We will employ the linear sigma model, the most popular one in 
this context, which seems to have the essential attributes of being simple but realistic 
enough: it has the correct chiral symmetry properties and also appropriately describes 
the low energy phenomenology of pions. We will choose the parameters of the model 
to give reasonable values for three experimentally determined quantities : the mass 
of the pion m^, the pion decay constant /„. and the s wave tt — tt phase shifts above 
threshold. 

We remark that the theory only makes sense at cutoffs below the Landau pole. 
This limits the size of the renormalized coupling constant 
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where A is the bare coupling and II is the polarization defined below. Since the mass 
difference between the a and tt is directly proportional to X r (q 2 = rn a 2) } this leads to 
an upper bound for the the mass of the a resonance as a function of A. Therefore, 
unlike at tree level, the mass of the a in the quantized theory is constrained in this 
model. A is also constrained from the physical considerations that we want the mass 
of the a to be less than the cutoff. However, the mass of the a resonance increases 
as we decrease the cutoff since the renormalized coupling increases. This pins down 
the cutoff to lie between 0.7GeV and lGeV. For a cutoff of lGeV we find that 
\™ ax (q 2 = 0) ~ 11.5 (taking into account finite cutoff effects). In our numerical 
simulations we use slightly smaller X r in order to avoid numerical problems that 
occur in the vicinity of the Landau pole. 

2. The linear cr-model in the large N approximation 

We will use the 0(4) linear sigma model to describe the evolution of the pions. 
We are well aware of the limitations of this approach, that provides a reasonable 
phenomenological model only for a limited range of energies (typically smaller than 
1 GeV). In spite of its shortcomings, this model captures some of the essential physics 
involved in the dynamics of the phase transition that may produce disoriented chiral 
condensates. In particular, the chiral phase transition takes place at a reasonable 
temperature (T c = a/3/tt ), and the low energy tt - tt scattering amplitudes are 
reasonable in this model. The mesons are organized in an 0(4) vector $ = (<r, tt) and 
the action is (in natural units % = c = 1) 
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where we have use the Bjorken and Drell metric: (1, —1, —1, —1). We will describe 
the evolution of the mean value $ =< $ > and the two point correlation functions 
including the effects of quantum and thermal fluctuations. The complexity of the 
problem forces us to adopt some approximations. Perturbation theory is useless for 



our purpose 14 > 15 5 and a scheme which is non-perturbative in A must be adopted. In 
this context, the most popular approaches which allow for a real time analysis are the 
the Hartree (or Gaussian) ansatz 16 and the large N expansion of the O(N) sigma 
model 17 . We will adopt the latter since it presents several advantages. On the one 
hand, the expansion is systematic and allows us to study higher order corrections 
(work is in progress in this direction and results will be presented elsewhere 18 ), and 
on the other hand, when using the Hartree ansatz in the context of the study of 
DCC's, one is forced to take the large N limit as well(see Ref. 12 ). This is due to 
the well known fact that the Gaussian approximation violates the Goldstone theorem 
giving an unphysical (and not necessarily small) mass to the pions in the H = 
limit. Of course, the approximation we adopt here is not expected to capture all the 
features of the phase transition (as is well known, mean field theory fails to predict 
the correct critical exponents, but allows us to explore the strong coupling regime). 

The large N effective equations can be obtained in a variety of ways, which are 
extensively discussed in the literature 17 . A very convenient method is to use an 
alternative equivalent action, which is a functional of the original fields $ and of the 
auxiliary constrained field x = A($ 2 /2iV — v 2 ) 19 . 

X ] = f d 4 x{-^( a + x)<f>i + ^ + \xv 2 + Ha}. (3) 

As this action is now quadratic in $ we can perform the functional integral over those 
fields, and are left with a functional integration over x which, to leading order in 1/N 
can be calculated by the stationary phase method. The effective action, which is a 
functional of the mean fields of $ and x, is obtained by Legandre transformation of 
the generating functional of the connected Green's functions. To the lowest order (in 
large N expansion) the effective action is given by 

r[$, x] = j d A x{- l -UU + x)^ + ^ + l - X v 2 + Ha+ l -Ntr lnG" 1 }, (4) 

where 

G?(x,y) = i[n + x (x)]6 4 (x-y), (5) 

and where for notational convenience we drop the overbars and denote the expectation 
values as $ and x- 

From this effective action we derive the equations of motion for the fields $ and 
the equation of constraint (gap equation) for the composite field x which plays the 
role of the effective mass for the mean values 

(n x + x (x))<S> l (x) = HS ll (6) 



x (x) = X(-v 2 + $ 2 (;c) + NG (x, x j). 
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The function Go(x,x) that appears in (7) is the coincidence limit of the propaga- 
tor Go(x,y) that inverts the operator Gq 1 defined in (5). We can use an auxiliary 
quantum field <f>(x) where < <f>(x) > = and 
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to construct this propagator as Go(x } y) =< T (f>(x)cf>(y) >. 

The initial value problem associated with equations (5)-(7) will be solved in the 
next section. Here, we would like to address the issue of how to use this model to 
make contact with the phenomenology we want to describe. Thus, we must fix the 
values of the parameters appearing in the above equations so that they describe low 
energy pion physics. The measurable quantities we want to reproduce are the pion 
mass m v = 135 MeV } the pion decay constant /„. = 92.5 MeV as well as the s 
wave, 1 = phase shifts in the energy range 300 — 420AfeV . To fit these physical 
quantities we analyze our equations in the "true vacuum" state (i.e., in equilibrium 
at zero temperature). In this state the derivatives of the expectation values vanish, 



and we have $ 



X = Xv where u v and \ v are some constants whose values 



will be determined below. The physical masses can be related to the parameters of 
the theory by computing the inverse propagators of the pion and sigma fields. The 
s-wave phase shifts are determined from the tt — tt scattering amplitudes obtained 
in this approximation which are controlled by the exchange of the composite field x 
propagator. In the vacuum sector (where x = = and a = f v ) the 7r, <t, and 

X propagators are obtained by inverting the matrix inverse propagator 
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ith the 5 dimensional vector ^ = ($, x), and are given by 
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We use the notation that [d n q] = d n q/(2Tr) n . (Notice that in the initial value problem 
one can have 7r 8 - ^ 0. Examining Eq. (10) we realize that 6i 1 Go(x } y) is the pion 
propagator only when < TTi(x) >= 0.) 



The pion-pion scattering invariant T-matrix is given by: 

T = 6 l3 6 kl A{s) + 6 lk 6 3l A{t) + SuS jk A(u), (14) 

where the isospin indices are coupled to the four momenta pi } p 2} p3 7 p4 such 

that s = (pi + P2) 2 , t = (pi — P3) 2 u = (pi — P4) 2 . In leading order in large N the 
amplitudes A(s), Ait) and A(u) are exactly the x propagator, namely 

A(s) = -G xx (p 2 = s). (15) 

The large N expansion preserves the current algebra so that the usual low energy 
theorem is exact. That is, for small s,m 2 one easily shows that 
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This amplitude is independent of the coupling constant A. We thus need to go above 
threshold to determine the coupling constant. The 1 = scattering amplitude is 

A = 3A(s) + A(t) + A(u). (16) 

The s-wave scattering amplitude is obtained by integrating the the 1=0 scattering 
amplitude over angles. 



f l=0 = e^sin 6(3) = -U/l - — i I" dz A , (17) 
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where z = cos 6 and 6 is the scattering angle in the s channel center of mass system. 
The vacuum expectation value of a is determined by /„. 

<Tv = U- (18) 

On the other hand, in the vacuum we have tt = 0, and the pion inverse propagator is 
Gjj „(i,i/) = Gq (x,y)Sij. Therefore, the vacuum expectation value of x i 



is 



Xv = m v = m . (19) 

The sigma mass can be approximately determined in terms of the inverse sigma 
propagator as the zero in the real part of the inverse propagator (a more precise 
determination which gives a slightly different result is from the peak in the 7 = 0, 
/ = scattering amplitude). This leads to the equation: 

< = m 2 + a 2 Re[G xx (m 2 a )}. (20) 

We do not expect the theory to be valid for energies above 1 GeV since in that 
regime the correct dynamics is described by QCD. In fact if we raise A above 1 GeV, 
the maximum value of the renormalized coupling goes down, and the a mass becomes 



unacceptably low. Reasonable values for the mass of the a constrain A to be in the 
range 0.7GeV < A < lGeV. In that range the best value of the renormalized coupling 
X r is between 7-10. For this range of values for X r this model agrees qualitatively with 
low energy scattering data as was mentioned before. With this choice of parameters, 
it is difficult to obtain values of the a mass higher than 450AfeV. 

3. Cooling by expansion 

To analyze the possibility of forming DCC's we should take into account the spe- 
cific features characterizing the situation after a highly energetic heavy ion collision. 
Experimentally, a flat plateau in the distribution of produced particles per unit ra- 
pidity is observed in the central rapidity region (these results are obtained in p-p 
and other collisions). This suggests the existence of an approximate Lorentz boost 
invariance. Thus, the simplest picture of a collision, due to Landau 25 , is one in which 
the excited nuclei are highly contracted pancakes receding away from the collision 
point at approximately the speed of light. The boost invariance implies that the 
evolution of the "hot plasma" that is left in between the nuclei looks the same when 
viewed from different inertial frames. Of course, this is an approximate picture that 
is not valid for large values of the spatial rapidity and for transverse distances of the 
order of the nucleus size. The existence of this approximate symmetry can be used 
to make a very simple hydrodynamical model 13 that, in some cases, may describe 
the evolution of the plasma. It is worth reviewing these ideas very briefly. Firstly, 
one should recognize that the natural coordinates to make a boost invariant model 
are the proper time r and the spatial rapidity 77, defined as 



The observed symmetry will be respected by the model if one imposes initial values 
on a r =constant hypersurface (and not at constant laboratory time t). 

According to this simple model of a collision, the plasma evolves in a highly 
inhomogeneous way when viewed from the laboratory frame. In fact, by analyzing 
a constant t surface we realize that the field configurations strongly depend on the 
spatial coordinate x. Near the light cone \x\ = t the system is "hot" (corresponding 
to small values of the proper time r). On the other hand, for small values of x (that 
correspond to larger values of r) the system is "colder". This type of configuration, 
hot in the outside and cold in the inside, is schematically known as "Baked Alaska". 

In this paper we want to study the formation of DCC's using some of the ideas 
presented above. We will not assume a quasi equilibrium situation or use a phe- 
nomenological hydrodynamical (or kinetic) model to describe the evolution of our 
system. On the contrary, we will study the nonequilibrium evolution in its full glory 
and solve equations (5)-(7), which include thermal and quantum effects. Using the co- 
ordinates (21) and fixing boost invariant initial conditions at an initial proper time r 
we introduce the expansion and "cooling" of the plasma in a natural way. Therefore, 
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we do not need to introduce any ad hoc cooling mechanisms by hand. The cooling, if 
any, will appear as a result of the evolution, which is fully out of equilibrium (a sim- 
ilar approach in different kind of expansion has been investigated independently by 
Gavin and Mueller 11 ). In this way, we can really test the validity of the "quenching" 
approximation that has been almost universally used when analyzing the evolution 
of DCC's. 

3.1. The equations 

In the coordinates (21) Minkowsky's arc element is 

ds 2 = dr 2 - r 2 dn 2 - dx\. (22) 

To study the evolution of the mean fields and correlations in these coordinates, our 
first task is to rewrite equations (5)-(7) using the new variables. Assuming that 
the mean values $ and x are functions of r only (homogeneity in the constant r 
hypersurface) we have 

r~ x d T rd T $,-(r) + x(r) $,-(r) = H6 lX (23) 

x (r) = \(-v 2 + $?(r) + N < <t> 2 (x, r) >), (24) 
where the quantum field ^(x,r) satisfies the Klein Gordon equation 

(r^drTdr -T- 2 d 2 -dl+ X (x))^X,T) = 0. (25) 

The quantum field <f>(x, r) defined here is an auxiliary field which allows us to calculate 
the Wightman function Go(x } y) by taking the expectation value: 

G (x, y; t) =< (f>(x, t) (f>(y, r) > . (26) 

When the expectation value of the pion field is zero, then this field corresponds to 
one component of the pion field. 

We expand this field in an orthonormal basis 

</>(??, £_l,t) = J[d 3 k](exp(ik ■ x)/ k (r) a k + h.c), (27) 

where k-x = k^n -\- k±x±, [d 3 k] = dk r] d 2 kj_/(2ir) 3 and the mode functions /^(t) evolve 
according to: 

A + (^ + kl + X(r) + t^)A = 0. (28) 



A dot here denotes the derivative with respect to the proper time r. The expectation 
value < <P 2 (x } t) > can be expressed in terms of the mode functions / k and of the 
distribution functions 

n k =< a[a h >, g h = < a k a k >, (29) 

which entirely characterize the initial state of the quantum field. For simplicity, we 
will assume that the initial state is described by a density matrix which is diagonal 
in the number basis (like a thermal state). In such a case, the only nonvanishing 
distribution is n k . Thus, replacing the above expressions in (24) we have: 

X (r) = A!;-*; 2 + $ 2 ( T ) + l -N J [d\]\f k (r)\ 2 (1 + 2 n k )). (30) 

We assume that the initial density matrix at r is one of local thermal equilibrium in 
the comoving frame 

n k (r ) ^ n k = ^1 _ x , 

where ft = 1/T , E° = ^fk 2 /r 2 + k\ + x (r ). 

To renormalize the theory we absorb the quadratic divergences in the bare mass 
Xv 2 and the logarithmic ones in the coupling constant A. After a few simple manip- 
ulations (that involve adding and subtracting the appropriate terms in (30)) we can 
write the equation for x as 

Xfj)=x(r ) + |(*?(r)-<I»?(T- )) + ^y[A]{|/ t (r)| 2 (l + 2, !t )- 
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where cj k (r) = (A; 2 /r 2 + k\ + x( T o)) 1/ ' 2 an d the renormalized coupling A r and Z are 
defined as 

N r A 1 

A = A r /Z, Z = l-X r 6X, 6X = — [d\]-T— . (33) 

4r J u£(t) 

The initial value x( T o) comes from solving the gap equation (30) at the initial time 
using the fact that in the initial equilibrium state at r , <& eq (T ) = H6n/x e q(To). 
Notice that the set of coupled equations, (23), (28), (32), are entirely written in 
terms of the renormalized quantities A r and xo- Changing the value of the cutoff will 
change both the value of the integrals and the value of Z in Eq. (32) leaving x cutoff 
independent. 
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3.2. The correlation function 



Expanding the field <f) in terms of the first order adiabatic mode functions, 
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where cu k (r) = + k\ + xM) 1/2 , we he 

<f>(ri, x u T ) = -L f [J 3 k](exp(zk • x)/ k °(r) a k (r) + h.c). (35) 



The two sets of creation and annihilation operators are connected by a Bogoliubov 
transformation: 

a k (r) = a(k,r)a k + /3(k,r)a[ k , (36) 
where the Bogoliubov coeeficients are given by 

a(k,r) = *(^^-^/ k ) (37) 

/3(k,r) = z( /k °^-^ /k ). (38) 
The time dependent creation and annihilation operators satisfy 

«k/ k ° + 4/k° = (39) 

in order for the usual canonical commutation relations to hold. Then we can write 
the Green function G° as follows: 

G°(x,y;r) = i/[J 3 k] e *-( x -y)|/ k (r)| 2 (l + 2n t (r )) 

= i /[j3 k] , k -(x-v) (l + 2 n k (r) + 2Re[ g ,(r)e- 21 ^ ^]) ^ 
t J 2cu k (r) ' 

where 

<4(r)a q (r)> = (27r) 3 «5 3 (k - q)n k (r) 
< a_ k {T) a(i {T) > = (27r) 3 «5 3 (k-q) 5k (r). (41) 

In terms of the initial distribution of particles nfc(r ) and the Bogoliubov coefficients 
a and f3 we have: 

n k (r) = n k (r ) + |/3(k, r)| 2 (l + 2n k (r )) (42) 



</ k (r) = a(k,r)/3(k,r)(l + 2n k (To)). 
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The variables n(r) and g(r) now have the physical interpretation as the interpolating 
phase space number and pair density for the case where < 7r 8 - >= 0. These operators 
coincide with the time independent number and pair densities defined by (29) at r , 
and in the out regime become the physically measurable number and pair densities. 
For < TTi >^ the actual pion two point function is more complicated than G° and 
this is only one piece of that Green's function. Of course in the vacuum Isospin 
conservation requires < 7r 8 - >= 0, so that once we are approaching the final true 
vacuum state during the cooling process, one can use these interpolating number and 
pair operators to describe the physics of the problem. Also note that during the 
evolution, uj\ can become negative (the instability phase) and for those modes the 
adiabatic basis does not exist. 

4. Initial conditions and results 

In order to solve equations (23), (28) and (30) as an initial value problem, we 
need to give Cauchy data (the function and its derivatives) for the mean values $;(t) 
and for the mode functions f\{r). It is possible to choose the initial data /k(^"o) and 
fk(To) so that the vacuum state coincides with the ordinary Minkowsky vacuum, at 
least for high momentum. It can be shown that this is accomplished by taking the 
"zero order adiabatic" vacuum 26 where 



We should point out that the initial value of \ wm always be restricted to being 
positive. As we said, x( T ) 1S the effective mass of the quasiparticles. Therefore, taking 
a negative initial value for \ implies that we are "turning the effective potential upside 
down." In our view, this cannot be the consequence of forming a "peculiar" initial 
state but must rather be the consequence of the cooling mechanism, which is entirely 
produced by the expansion. In fact, starting with a negative Xo 1S what is done when 
studying this problem by using the quench approximation: one starts with a hot 
initial state and lets it evolve in the low temperature effective potential. It should 
be clear by now that this is drastically different from our approach. We will study 
the selfconsistent evolution of \ starting from a "hot" initial value and follow its 
evolution. 

We want to start in the quark - gluon plasma phase above the chiral phase 
transition which places constraints on the initial energy density present at r . We 
then assume that slightly above this transition it is reasonable to model the chi- 
ral transition with the effective Lagrangian of the sigma model. Not knowing ex- 
actly what this proper time is, we will consider here reasonable initial proper times 
(lfm/c < r < 4/ra/c) and study the effect of the initial proper time r on the 
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where cu k (r) = (PJt 2 + k\ + x{r)) 1/2 - 



production of instabilities. Taking larger proper times as the starting point for our 
calculation reduces even further the possibility of instability growth. 

The first issue we will examine here is the existence of instabilities. Notice that 
there are unstable modes fk whenever (& 2 + l/4)/r 2 + k\ + x is negative so that only 
the long wave length modes can become unstable. x can classically become negative 
only when 



so that the most unstable case has $(r) = 0. Therefore, at the classical level the bare 
coupling must satisfy 



for there to be any instabilities. Once instabilities grow, they cause exponential 
growth of the modes in the mode sum which contributes a positive quantity to the 
equation for x- Thus, the stronger the coupling the quicker x returns to a positive 
value. 

In previous works the presence of unstable modes was assumed by imposing the 
quench approximation. In our numerical studies we find that, due to the strong 
coupling, in a wide class of initial conditions consistent with reasonable fluctuations 
found in a thermal distribution no instabilities develop. In our simulations the one 
thing we did not change was the value of the composite field x which was fixed to be 
the solution of the gap equation in the initial thermal state. We also maintained the 
constraint that the initial values of the expectation values of the field, namely 7t(t ) 
and cr(r ) satisfied the relationship: 



where is the equilibrium value of $ at the initial temperature T . This last con- 
straint is reasonable since it doe not cause much energy to make a rotation in the 
direction of 7r 4 -. Thus we chose different initial a and 7r 8 - consistent with the above 
constraints and then probed the effect of different initial values for a and 7r 4 -. Since the 
results for 7r 8 - ^ were similar to those when (7^0 we mainly present here the results 
for the latter case. We have also surveyed other possibilities that violate the above 
constraints such as choosing initial a = 0. In this case we found that instabilities 
(even with a ^ 0) have shorter life times than those with the above constraint. 

We have performed numerical simulations on the connection machine CM-5 using 
a grid that has 100 modes in the transverse direction and 100(t/t o ) modes in the r\ 
direction with dk^/fm = A/100, dk ± = A/100. 

First, let us consider the case where the initial value cr(r ) is the thermal equi- 
librium one corresponding to a temperature of 200MeV . We varied the value of the 
initial proper time derivative of the sigma field expectation value and found that there 
is a narrow range of initial values that lead to the growth of instabilities, namely 
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0.25/m" 2 < \&\ < 1.3/m" 2 . 



Surprisingly when \a\ > 1.3/m~ 2 instabilities no longer occur. 

Figures 1-2 summarize the results of the numerical simulation for the evolution of 
the system (23), (28) and (32). We display the auxiliary field x m units of fm~ 2 , the 
classical fields $ in units of jm~ x and the proper time in units of jm~ x . In Fig. 1 and 
Fig. 2 the proper time evolutions of the auxiliary field x field are presented for two 
different values of /„., where the initial conditions were fixed at r = 1. In the case of 

= 92.5MeV we see that the instability lasts for less than 3/ra. As discussed earlier, 
the regime of exponential growth of the unstable modes occurs whenever x < 0. We 
notice that when we rotate the expectation value from the a to the tt 1 direction 
initially, the typical time that instability exists is of the same magnitude. We also 
notice that if we choose the time derivative to be zero and start from an equilibrium 
configuration, then the expansion alone is insufficient to generate instabilities. We 
find that in order to generate instabilities we require fluctuations in the classical 
kinetic terms such as a or 7r 4 -, and as discussed earlier, these initial conditions must 
be in a very narrow range to produce instabilities. Thus the rapid quench conditions 
assumed by other authors comprises only a small region of the phase space of initial 
fluctuations expected in an initial thermal distribution. As expected, at late proper 
times, the auxiliary field approaches its equilibrium value of m\. For the case where 
/,r = 125MeV } a value favored by fitting the low energy scattering, we see a very 
similar pattern, showing that our main results concerning the small range of initial 
conditions that lead to instabilities are not affected by a 30% change in the value of 
this parameter. 

A crude estimate for the size of a disoriented chiral condensate is to multiply a 
typical life time of the instabilities by the speed of light. With this estimate the 
size of these regions are of the order of a few fermi. Of course one needs to study 
the growth of inhomogeneous instabilities to make any definite statements about this 
size. We were hoping that the correlation functions we have calculated could be 
interpreted easily in terms of a length parameter associated with the size of these 
regions. However, as we will see below, this didn't occur. In Fig. 3 we plot the 
proper time evolution of the classical fields a and tt for the same two choices for /„. as 
in Fig. 1 and Fig. 2. We see that the sigma field asymptotes to its vacuum value /„. 
and the tt field gradually converges to its equilibrium value of zero. In Fig. 4 we show 
the effect of starting the initial value problem at later times, namely r = 2.5, 4 /ra, 
and compare them to the case previously studied with r = which had a modest 

region of instability growth. We see that as we increase r we decrease the possibility 
of instability growth, and with these late initial times it is not possible to produce 
instabilities even with kinetic energy fluctuations. In Fig. 5 we study the effect of 
the initial temperature on our time evolution problem. For /„. = 92.5 the critical 
temperature is 160MeV (in the absence of explicit symmetry breaking term). We 
see that in the vicinity of the critical temperature the effects of varying the initial 
temperature is minor. In Fig. 6 we study the proper time evolution of the effective 
number density for various initial conditions. In thermal equilibrium one expects for 
an isentropic expansion in boost invariant coordinates that st = constant. Since the 
number density is proportional to the entropy density one expects that once particle 
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Figure 1: Proper time evolution of the x field for four different initial conditions with 
f n = 92.5MeV. 



production stops that n(r)r — > constant. We see this trend in this figure, showing 
that we are reaching the out regime as the system expands. The breaks in some 
of these curves at early values of r are a result of the fact that the interpolating 
number density cannot be defined when \ 1S negative. Next, we are interested in 
seeing how our results differ from the case where the system evolves in local thermal 
equilibrium. This is described by two correlation lengths, the inverse of the effective 
pion mass associated with x( T ) ■> an d the inverse of the proper time evolving effective 
temperature T(t) = T(t )(t /t) 3 corresponding to a boost invariant hydrodynamical 
models. We see from Fig. 7 that in the case of cr(l) = <tx, 7r*(l) = and <r(l) = — 1, 
where maximum instability exists, complex structures are formed as contrasted to 
the local thermal equilibrium evolution. The interpolating phase space distribution 
n(^,kj_,r) obtained numerically, clearly exhibits a larger correlation length in the 
transverse direction than the nx(^,kj_,r) equilibrium one, and has correlation in 
rapidity of the order of 1-2 units of rapidity. We notice that in both directions there 
are structures that does not lend themselves to a simple interpretation. On the 
other hand, the local thermal equilibrium evolution is quite regular apart from the 
normalization of the distributions that are changing with time due to oscillation in 
the quantity x( T ) which is damped to its equilibrium value once the system expands 
sufficiently. In Fig. 8 we look at the time evolution of the interpolating number 
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Figure 2: Same as Fig. 1, but for /„. = 125MeV. 



operator for the case in which we did not have any fluctuations in the kinetic energy. 
Here we are still far from equilibrium, and we see that unlike the equilibrium case the 
correlation length in rapidity space is not decreasing with proper time. In this case 
where there are no instabilities we do not see complicated structures. The transverse 
distribution is similar to the equilibrium case. The pair density function kj_, t) 
is even more elusive to parametrize than the single particle distribution function. In 
the case where there are no instabilities (<r(l) = 0) we see in Fig. 9 that although the 
transverse distribution is relatively simple, the distribution in the rj direction (whose 
fourier transform gives the rapidity distribution) has many length scales. When we 
have instabilities (<r(l) = —1) then both distributions are complicated and possess 
several length scales, as seen in the lower part of Fig. 9. We also directly computed 
the fourier transform of the pion correlator (the hrst line of Eq.(40) for the case of 
TTi = 0, and as expected, it has even more complex structure than the distributions n 
and g. In the case of non zero 7r 8 - fields we expect that the correlator will have similar 
structure to the case described above, and it will involve at least 3 time evolving length 
scales related to the inverse mass of the pion, the inverse of the effective temperature 
and possible length scales describing domain growth. 




Figure 3: Proper time evolution of the a and tt fields for the following initial condi- 
tions: Solid line is for ail ) = <tx, 7r* (1) = and <r(l) = f . Dotted line is for ail ) = <tx, 
7T* (1) = and <t(1) = 0. Dashed- dotted line is for cr(f) = <rx, 7r* ( 1 ) = and 
<j(l) = -1. Dashed line is for <r(l) = 0, tt*"(1) = °"x and <r(l) = -f . At T = 200MeV, 
a T = 0.3/m" 1 for f n = 92.5MeV and a T = 0.5/m" 1 for f n = 92.5MeV. 




Figure 5: Proper time evolution of the x field for three different initial thermal 
distributions with T = 200, 164, 150MeV for the case of <r(l) = <j t , tt 8 (1) = and 
a(l) = -l. 




Figure 6: Proper time evolution of the produced particle density nr = dN/dr]dx± for 
the same evolution shown in Fig. 1. 
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Figure 7: Slices of k v = and p = |kj_| = of the proper time evolution of the 
interpolating phase space particle number density n(k V7 kj_, r) for cr(l) = <rx, 7r* ( 1 ) = 
and <t(1) = — 1 compared with the corresponding local thermal equilibrium densities 
n T (^,k_L,r). 
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Figure 9: Slices of k v = and p = |kj_| = of the proper time evolution of the 
real part of the phase space pair density kj_, t) for cr(l) = <rx, 7r* ( 1 ) = and 

<t(1) = 0,-1. 



5. Conclusions 

In this paper we have performed numerical simulations in the regime of the chiral 
phase transition in the linear sigma model for a wide variety of initial conditions 
starting above the critical temperature for an expanding plasma of pions and sigmas. 
Assuming that this model gives a reasonable description in this temperature range, 
we found initial conditions where instabilities grew in the scenario required for the 
formation of disoriented chiral condensates. The constraints on our model which are 
imposed in order to approximate the low energy physics of pion interactions require 
a large renormalized coupling constant. This had the effect of rapidly damping the 
instabilities, and we found no evidence in this model for large domains of disoriented 
chiral condensates. We did, however, see rather large departures in the phase space 
number density from one which would result from an evolution in local thermal equi- 
librium. These departures show a narrowing of the momentum space distribution 
which could be interpreted as a larger spatial correlation length. However we did 
not find a simple way of extracting correlation lengths from our results. In our sim- 
ulations we assumed a reasonable mechanism for cooling- namely the expansion of 
the plasma following its production in a collision. In the future we hope to include 
scattering effects which will introduce another time scale, namely the equilibration 
time scale, into the problem, as well as study the nonlinear sigma model to see if the 
results differ significantly from those found here. 
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